Spherical excision for moving black holes and 
summation by parts for axisymmetric systems 
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It is expected that the realization of a convergent and long-term stable numerical code for the 
simulation of a black hole inspiral collision will depend greatly upon the construction of stable 
algorithms capable of handling smooth and, most likely, time dependent boundaries. After deriv- 
ing single grid, energy conserving discretizations for axisymmetric systems containing the axis of 
symmetry, we present a new excision method for moving black holes using multiple overlapping 
coordinate patches, such that each boundary is fixed with respect to at least one coordinate system. 
This multiple coordinate structure eliminates all need for extrapolation, a commonly used procedure 
for moving boundaries in numerical relativity. 

We demonstrate this excision method by evolving a massless Klein-Gordon scalar field around 
a boosted Schwarzschild black hole in axisymmetry. The excision boundary is defined by a spher- 
ical coordinate system co-moving with the black hole. Our numerical experiments indicate that 
arbitrarily high boost velocities can be used without observing any sign of instability. 
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I. INTRODUCTION 

Inspiraling black holes arc among the strongest astro- 
physical sources of gravitational radiation. The expecta- 
tion that such systems may soon be studied with gravi- 
tational wave detectors has focused attention on solving 
Einstein's equations for predictions of gravitational wave 
content. Although the Einstein equations present sev- 
eral unique challenges to the numerical relativist on 
several of which we do not elaborate here, black holes 
present one particular additional challenge: they contain 
physical curvature singularities. While the infinities of 
the gravitational fields associated with this singularity 
cannot be represented directly on a computer, the space- 
time near the black hole must be given adequately to 
preserve the proper physics. 

Different strategies have thus been developed to com- 
putationally represent black holes, while removing the 
singularity from the grid. One method exploits the gauge 
freedom of general relativity by choosing a time coor- 
dinate that advances normally far from a singularity, 
slows down as a singularity is approached, and freezes 
in the immediate vicinity. Coordinates with this prop- 
erty are "singularity avoiding" While singularity 
avoiding coordinates have some advantages, one poten- 
tial disadvantage is that the hypersurfaces of constant 
time may become highly distorted, leading to large gra- 
dients in the metric components. These slice-stretching 
(or "grid-stretching") effects, however, can be partially 
avoided through an advantageous combination of lapse 
and shift conditions. For example, long-term evolutions 
of single black holes have been reported by Alcubierre 
et al. |5(. Singularity avoiding slicings may be combined 
with black hole excision, a second method for removing 
the singularities from the computational domain. Cur- 
rently, long-term binary black hole evolutions have only 
been performed using both techniques together. 

Excision is based on the physical properties of event 



horizons and the expectation that singularities always 
form within such horizons, and thus cannot be seen by 
distant observers, as formulated by the Cosmic Censor 
Conjecture 0. As no future-directed causal curve con- 
nects events inside the black hole to events outside, Un- 
ruh proposed that one could simply remove the black 
hole from the computational domain, leaving the exte- 
rior computation unaffected 0. Thus the black hole 
singularity is removed by placing an inner boundary 
on the computational domain at or within the event 
horizon. Excision has been extensively used in numer- 
ical relativity in the context of Cauchy formulations 
mmSIIllIlllllIllIllIllIllIi In particular, 
excision with moving boundaries, which is the primary 
focus of this paper, was explored in [T^ Hsl IT^ . 

The physical principles that form the basis of exci- 
sion make the idea beautiful in its simplicity. Translat- 
ing them into a workable numerical recipe for black hole 
evolutions, on the other hand, requires some attention to 
detail. Two general questions arise regarding the imple- 
mentation of excision, (1) Where and how to define the 
inner boundary? and (2) How to move the boundary? 
The first question applies to all excision algorithms, while 
the last question is specific to implementations where the 
excision boundary moves with respect to the grid. In ad- 
dressing these questions we assume a symmetric (or at 
least strongly) hyperbolic formulation [l^. This is be- 
cause excision fundamentally relies on the characteristic 
structure of the Einstein equations near event horizons, a 
structure which can only be completely defined and un- 
derstood for strongly and symmetric hyperbolic sets of 
equations. 

The first question involves several considerations, in- 
cluding the location of the boundary, its geometry, and its 
discrete representation. The requirement that all modes 
at the excision boundary are leaving the computational 
domain can be non trivial. It may appear that this condi- 
tion would be satisfied simply by choosing any boundary 
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within the event horizon (or, for practical purposes, the 
apparent horizon). However, the outflow property of the 
excision boundary depends on the characteristic speeds 
of the system in the normal directions to the bound- 
ary. For example, in the analytic Schwarzschild solu- 
tion, assuming that the system has characteristic speeds 
bounded by the light cone, a spherical boundary can be 
excised at r < 2AI. A cubical boundary, on the other 
hand, imposes an onerous restriction on the excision vol- 
ume: in Cartesian Kerr-Schild coordinates the faces of 
a cube centered on the black hole must be less than 
4V3/9Af « 0.7698Af in lengthfll l^^. Remarkably, 
as was first noticed by Lehner |2j , for the rotating Kerr 
solution in Kerr-Schild coordinates a well-defined cubi- 
cal excision boundary is impossible for interesting values 
of the spin parameter. (See the appendix for further 
discussion.) Whereas with a pscudospcctral collocation 
method the implementation of a smooth spherical ex- 
cision boundary is trivial |22| . this is generally not the 
case for finite differencing. As may be expected, smooth 
boundaries, which can be adapted to the spacetime ge- 
ometry, allow the excision boundary to be as far from the 
singularity as possible, making the most efficient use of 
the technique. 

The discrete representation of boundaries can be a 
delicate issue, especially in numerical relativity where 
many large-scale finite difference computations are done 
in Cartesian coordinates. We focus our attention on 
smooth boundaries that may be defined as a constant 
value in the computational coordinates, e.g., r = rg in 
spherical coordinates describes a simple spherical bound- 
ary. The importance of accurately representing smooth 
boundaries has been demonstrated for the Euler equa- 
tions, for example, by Dadonc and Grossman for finite 
volume methods , and Bassi and Rebay using fi- 
nite elements. Bassi and Rebay studied high resolution 
planar fiuid fiow around a cylinder. They report spuri- 
ous entropy production near the cylinder wall, which cor- 
rupts the solution even on extremely refined grids, when 
the cylindrical boundary is approximated by a polygon. 
Furthermore, in the conformal field equations approach 
to general relativity, a smooth boundary is required to 
avoid uncontrollable numerical constraint violation |25j . 

The second question applies to excision boundaries 
that move with respect to the grid. When the inner 
boundary moves, points that previously were excised en- 
ter the physical part of the grid, and must be provided 
with physical data for all fields. In recently proposed 
excision algorithms, these data arc obtained by extrapo- 
lating the solution from the physical domain of the cal- 
culation. Numerical experiments have indicated that the 
stability of the method is very sensitive to the details of 
the extrapolation, see e.g., Ref. 0,0,11^. To exam- 
ine the black hole excision problem with moving inner 
boundaries, we adopt an approach with some unique fea- 
tures. The heart of our method for moving excision is to 
use multiple coordinate patches such that each boundary 
is at a fixed location in one coordinate system. Adapting 



coordinate patches to the boundary geometry allows us 
to excise as far from the singularity as possible and sim- 
plifies the determination of the outfiow character of the 
excision boundary. The motion of the boundaries is in- 
corporated through the relationships among the various 
coordinate systems. The grids representing the different 
coordinate patches overlap and communicate via inter- 
polation. This technique is an extension of the one used 
in well-posedness proofs for problems in general domains 
(see Sec. 13.4 of |23|). In this paper we demonstrate the 
algorithm by solving the massless Klein-Gordon equation 
on a fixed, boosted Schwarzschild background. We find 
that the algorithm is stable for (apparently) all values of 
the boost parameter, /3 — v/c, and present results here 
showing stable evolutions for several cases with /? < 0.95. 

We specialize to axially symmetric spacetimcs to re- 
duce the computational requirements for our single- 
processor code. Axially symmetric spacetimcs have 
sometimes been avoided in numerical relativity, with no- 
table exceptions, see e.g., owing to the difficulties in 
developing stable finite difference equations containing 
the axis of symmetry. In this paper we further present 
finite differencing methods for the wave equation in axi- 
ally symmetric spacetimes in canonical cylindrical and 
spherical coordinates. These differencing schemes are 
second order accurate and their stability for a single grid 
is proved using the energy method |27l | . Maximally dissi- 
pative boundary conditions are applied using the projec- 
tion method of Olsson |2^. We present the differencing 
algorithm in detail, and indicate precisely how boundary 
conditions are applied. 

This paper is organized as follows: In Sec. ^] we mo- 
tivate our approach and review the overlapping grid 
method. We recall the concept of conserved energy 
for a first order symmetrizable hyperbolic system in 
Sec, mil and provide an energy preserving discretization. 
In Sec. IIVI we analyze the axisymmetric wave equation 
around a Minkowski background as an introduction to 
our numerical methods. The analysis is then repeated 
for the black hole background case in Sec. The ex- 
cision of a boosted black hole with the overlapping grid 
method is described in Sec. IVII The numerical experi- 
ments, along with several convergence tests, are included 
in Sec. Em 



II. OVERVIEW 

Our primary goal is to obtain a numerical algorithm 
for excision with moving black holes that is stable and 
convergent (in the limit that the mesh spacing goes to 
zero). These desired properties for the discrete system 
closely mirror the continuum properties of well-posed ini- 
tial boundary value problems (IB VPs): the existence of 
a unique solution that depends continuously on the ini- 
tial and boundary data. Furthermore, we believe that 
we will not obtain long-term convergent discrete solu- 
tions unless the underlying continuum problem is also 
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FIG. 1: A singularity S surrounded by an event horizon dO. is 
moving with respect to the base coordinate system. A coordi- 
nate patch (shaded region) adapted to 9f2 follows the motion 
of the singularity. With respect to this patch, dO. is a purely 
outflow boundary and requires no boundary conditions. The 
base system terminates somewhere inside the shaded region 
and it gets boundary data from the moving patch. Similarly, 
the data at the outer boundary of the moving patch is taken 
from the base system. 



well-posed. Unfortunately there are few mathematical 
results concerning the well-posedness of general classes 
of equations. The energy method, however, can be used 
with symmetric hyperbolic IBVPs, and gives sufficient 
conditions for wcU-poscdness. 

When a black hole moves with respect to some co- 
ordinate system, the inner excision boundary must also 
move. We use multiple coordinate patches, such that 
every boundary is fixed with respect to at least one co- 
ordinate system. Coordinate transformations relate the 
coordinate systems, and become time dependent when 
the hole moves. The movement of the inner boundary is 
also expressed by these time-dependent coordinate trans- 
formations. These ideas are illustrated in Fig. In 
our axially symmetric model problem of a scalar field 
on a boosted Schwarzschild spacetime, the computational 
frame is covered with cylindrical coordinates, while a sec- 
ond patch of spherical coordinates is co-moving with the 
black hole. (In these coordinates the event horizon is 
always located at r = 2M while the time coordinate is 
taken from the cylindrical patch so that data on all grids 
are simultaneous.) The inner boundary of the spher- 
ical grid, located at or within the event horizon, is a 
simple outfiow boundary, and requires no boundary con- 
dition. The cylindrical domain has an inner boundary 
somewhere near the black hole, whether inside or out- 



side of the horizon is immaterial, as long as it is covered 
by the spherical coordinate patch. An exchange of infor- 
mation between the two coordinate patches is required 
to provide boundary conditions at the inner cylindrical 
boundary and the outer spherical boundary. 

On each grid the discrete system is constructed using 
the energy method [23|- We define an energy for the 
semi-discrete system and, using difference operators that 
satisfy summation by parts, we obtain a discrete energy 
estimate Well-posed boundary conditions can then 
be identified by controlling the boundary terms of the 
discrete energy estimate. The conditions are discretized 
using Olsson's projection method In particular, the 
symmetry axis (p = in canonical cylindrical coordi- 
nates) is included in the discrete energy estimate 
allowing us to naturally obtain a stable discretization for 
axisymmetric systems. 

We implement our excision algorithm using overlap- 
ping grids, also known as composite mesh difference 
method [23, |^ 1^3 • The two grids are coupled by in- 
terpolation, which is done for all the components of the 
fields being evolved. If the system is hyperbolic this 
means that one is actually over specifying the problem. 
However, as it is pointed out in Sec. 13.4 of |23| and as 
it is confirmed by our experiments, this does not lead 
to a numerical instability. The fully discretized system 
is completed by integrating the semi-discrete equations 
with an appropriate method for ODEs; we choose third 
and fourth order Runge-Kutta, which does not spoil the 
energy estimate of the semi-discrete system |0 . Kreiss- 
Oliger dissipation is added to the scheme, as some ex- 
plicit dissipation is generally necessary for stability with 
overlapping grids |35j- Whereas the stability theory for 
overlapping grids for elliptic problems is well developed, 
there are very few results concerning hyperbolic systems. 
Starius presents a stability proof for overlapping grids in 
one dimension |3ll | . Finally, we note that Thornburg has 
also explored multiple grids in the context of numerical 
relativity with black hole excision [^Bl • 

The structure of the overlapping grids used in this work 
is illustrated in Fig.|Sl The additional complication of the 
axis of symmetry is discussed below. For simplicity we 
choose the outer boundary to be of rectangular shape. 
The introduction of a smooth spherical outer boundary, 
along with another grid overlapping with the base cylin- 
drical grid, is certainly possible and, we believe, likely to 
improve the absorbing character of the outer boundary 
when the incoming fields are set to zero. 



III. THE WAVE EQUATION 

To demonstrate our excision algorithm, we choose the 
evolution of a massless Klein-Gordon scalar field on an 
axisymmetric, boosted Schwarzschild background as a 
model problem. In this section we summarize basic defi- 
nitions for linear, first order hyperbolic initial-boundary 
value problems [l^ |23 . We employ the energy method 
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to identify wcU-poscd boundary conditions. The discrete 
version of this method, based on difference operators sat- 
isfying the summation by parts rule |37| . is then used 
to discretize the right hand side of the system and the 
boundary conditions on a single rectangular grid. (For an 
introduction to these methods in the context of numer- 
ical relativity sec Refs. [Ullllll.) We then introduce 
the axisymmetric scalar field equations on a curved back- 
ground, along with their semi-discrete approximation. 

In this paper we adopt the Einstein summation con- 
vention and geometrized units (G = c = 1). Latin indices 
range over the spatial dimensions, and Greek indices la- 
bel spacetime components. 



A. Hyperbolic systems in first order form 

Consider a linear, first order, hyperbolic IBVP in two 
spatial dimensions, consisting of a system of partial dif- 
ferential evolution equations, and initial and boundary 
data, of the form 

dtu = A'{t,x)diU + B{t,x)u {t,x) e [0,T] x 
u(0,f) = /(f) xen (2) 
Lu{t, x) = g{t, x) {t, x) e [0, T] x dil, (3) 

where i = 1,2, u = u{t,x) and f{x) are vector valued 
functions with m components, A'' and B arc m x m ma- 
trices that depend on the spacetime coordinates but not 
on the solution u, and di stands for d/d^i. The boundary 
of C is assumed to be a simple smooth curve. The 
operator L and the data g that appear in the boundary 
condition Q will be defined below in Eqs. (|10llll|) . 

1. Strong and symmetric hyperbolicity 

System (0-® is said to be strongly hyperbolic in O C 
[0,T] X il if, at each point {to,XQ) £ O, the matrix 



(4) 



uJi + UJ2 = 1, can brought into 



P{to,xa,L^) = A^{to,xo)i^j 

with w e and \uj\'^ 
real diagonal form by a transformation T{uj), such that 
T{u}) and T~^{uj) are uniformly bounded with respect to 
uj. The system is said to be symmetric or symmetrizable 
hyperbolic in O if, at each point {Lq, xq) G O, there exists 
a smooth, symmetric positive definite matrix H{tQ,XQ), 
independent of uj, such that HA^ — {HA' Y' for i — 1,2. 
The matrix H is usually called the symmetrizer. Clearly, 
a symmetric hyperbolic system is also strongly hyper- 
bolic. Strong hyperbolicity is a necessary condition for 
well-posedness and consequently for the construction of 
stable numerical schemes. 



2. Characteristic speeds 

The characteristic speeds in the direction ft = 
(71,1,712) G M^, with nl + n2 = 1, at the point (^0,^0) G 



[0, T] X J7 are the eigenvalues of A"(io, ^o) = f^i^H^o, xq). 
In Sec. IVIll we will show how the maximum value of the 
characteristic speeds in the region [0, T] x il can be used 
to compute an upper bound for the ratio between the 
time step and the spatial mesh size. 



3. Energy method 

The specification of proper boundary conditions re- 
quires careful consideration in order to achieve a well- 
posed IBVP, and we use the energy method to identify 
appropriate boundary conditions [2^ |40| . Here one de- 
fines the energy of the system at time t to be 

E{t) = \\u{t,-)\\% = / u^{t,x)H{t,x)u{t,x)(fx, (5) 
Jn 

where H is some positive definite m x m symmetric ma- 
trix and denotes the transpose of u. To ensure contin- 
uous dependence of the solution on the initial and bound- 
ary data, the energy must be bounded in terms of appro- 
priate norms of the data. To determine this bound one 
usually takes a time derivative of the energy (0), with 
the further assumptions that m is a smooth solution of 
^ and that iJ is a symmetrizer. The energy estimate is 
then 

^E{t) = / u^HA^'uda (6) 
dt Jon 

+ [ {dtH + HB + {HBf - d^{HA')) ud^x , 
Jn 

where Gauss' theorem was used to obtain the boundary 
term, and rij is the outward unit normal to the boundary 
To control the growth of the energy of the solution, 
we naturally need to control both the boundary and vol- 
ume integrals. We consider the boundary integral first. 

The matrix HA" is symmetric, and can be brought 
into diagonal form by an orthogonal transformation 
Qin), 

Q^{n)HA"Q{n) = A = diag(A+, -A„, 0) , (7) 

where K± > are positive definite diagonal matrices, 
the eigenvalues of which, in general, do not coincide with 
the characteristic speeds. This allows one to rewrite the 
integrand of the boundary integral in Eq. JH)) by intro- 
ducing the vector = (i„(+A+;")^u'(-^-'"), = 
Q^{n)u as the difference between two non- negative 
terms. 



u^HA^u 



^(-hA+;«)^^^^,(+A+;„) 



(8) 



The components of w^"'^ are the characteristic variables 
in the direction n. In particular, the components of 



w 



(+A+;n) 



are the ingoing characteristic variables, and 
the components of w^~^~'"^ are the outgoing charac- 
teristic variables. We see that prescribing homogeneous 
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boundary conditions (u'^"'"^+'"^ = S'w^"^"'"-', with S suf- 
ficiently small, i.e., S'^A+S < A_), ensures that the 
boundary term will give a non-positive contribution to 
the energy estimate. The S = case (no coupling) is of 
particular interest as it usually yields a good approxima- 
tion for absorbing (Sommerfeld) boundary conditions. 

The second term of the energy estimate the vol- 
ume integral, can be estimated by 2a||u(t, •)|||f , where 
a= imax(t,s) \\dtH + H B + {H B)^ - d,{H A')\\ is a con- 
stant that does not depend on the solution. Thus, for 
homogeneous boundary conditions we have 

jJuit,-)rH<2a\\uit,-)\\l, (9) 

which impUes that ||w(t, •)||i/ < exp(a<)||/||i/. Simi- 
lar energy estimates can be obtained for inhomogeneous 
boundary conditions 113,^3; i-e-, 

^(+A+;„) ^5^(-A_;„)^ ^^^ (10) 

where g has to satisfy compatibility conditions with the 
initial data. 

Boundary conditions of the form H10() are referred to 
as maximally dissipative boundary conditions |4ll | . From 
Eq. I|l()() we see that the operator L introduced in iPJ has 
the form 

L = pWQ^(n)-5P(-)g^(n), (11) 

where P^+'> (w^+\w^-\w°)'^ = (u;(+),0,0)^ and 
p(-)(u;(+)^iy(-)^u;0)T ^ (Q ^ w'^~\ O f . Finally, it is im- 
portant to recognize that if the matrix H is not a sym- 
mctrizer, as would be the case if H symmetrizes but 
fails to be positive definite, then the boundary condition 
above will not, in general, lead to a well-posed problem, 
as one would likely end up specifying boundary data to 
the wrong quantities. 

4- Strict stability 

Discretizing the spatial derivatives in the right hand 
side of system (J^l, but leaving time continuous, leads to 
a system of ODEs called the semi-discrete system. If an 
initial value problem satisfies the estimate ||u(t, •)||// < 
K e:>q>{at)\\u{0, ■)\\h at the continuum, it would be desir- 
able to obtain a discretization such that a similar esti- 
mate holds at the discrete level. Following |23, we will 
say that a semi-discrete system is strictly stable if 

Mm,, < Kse"''\\uiO)U (12) 

where as < a + 0{h) and || • \\h is a discrete energy 
consistent with the one of the continuum. 



5. Conserved energy 

Clearly, the requirement that HA^ be symmetric, does 
not uniquely determine the symmctrizcr. For example. 



if is a symmetrizer, then fH with / > is also a 
symmetrizer. In some circumstances, as for the scalar 
field considered here, it is possible to select a preferred 
symmetrizer which satisfies the additional requirement 

dtH + HB + {HBf ~di{HA')=^0. (13) 

When this condition holds, the energy defined by that 
symmetrizer will be conserved. By this we mean that 
the change in energy of our system is solely due to the 
boundary term of JB)), which can be controlled by using 
maximal dissipative boundary conditions H1U|) . In partic- 
ular, when homogeneous boundary conditions are used, 
or when no boundaries are present, the energy cannot 
increase. 

6. Energy conserving schemes 

Let us assume momentarily that that there exists a 
symmetrizer for which Eq. (|13|) holds (the system ad- 
mits a conserved energy), and that no boundaries are 
present. In the variable coefficient case (more precisely, 
if dj{HA^) 7^ for i = j), the naive discretization 
dfU = A^DiU -\- Bu, where u now represents a vector 
valued grid function, although strictly stable when a sec- 
ond order accurate centered difference operator is used 
[2^ , does not conserve the discrete energy 

E = {u, Hu)h = hih2 ^ ufjHijU.ij , (14) 

where Hij = H{t,Xij). Its time derivative gives 

^E[t) = {u, [HA\D,]u)h + (u, d,[HA')u)h + , (15) 
at 

where we used the fact that {v,Diu)h + {DiV,u)h — 0. 
The lack of a Leibniz rule at the discrete level is only 
partly responsible for this. In general, even if this rule 
were satisfied, the discrete estimate would not vanish, 
f^E = {u,{d,{HA') ~ D,{HA'))u)h ^ 0. Any semi- 
discrete approximation that preserves the discrete energy 
(|14|l is an energy conserving scheme. Remarkably, when- 
ever a system admits a conserved energy at the contin- 
uum it is always possible to construct an energy con- 
serving scheme [ll HI Is^. The following "1/2 + 1/2" 
splitting, for example, 

dtu = ^A'DiU + ^H-^D,{HA'u) (16) 

+ (^B-^H-'d,{HA')^u, 

ensures that the discrete energy H14|) remains constant. 
Clearly, an energy conserving scheme is strictly stable, 
since a = as = 0- We note that, depending on the 
problem, there may be alternative, simpler discretiza- 
tions than the "1/2 -I- 1/2" splitting which lead to the 
same energy estimate. Moreover, a discretization such as 
(I16|l is a consistent approximation of dtu = A^diU -t- Bu 
whether or not condition (|13|l holds. 
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7. Rectangular grid 

Consider a rectangular domain fl = {{x^ , x'^)\xl^-^^-^ < 
< ^Li^Ki^min < a:^ < xl^^J, with the grid points 

% = (^^min + 'i-'T'l^^min + J^2), « = 0, . . . , iVi and j = 

0, . . . , 7V2, and hk = (a;^^^ - xf„ij,)/A^fc, fc 1, 2. From 
the continuum analysis we expect that boundary data 
should be given to the incoming characteristic variables 
in the direction orthogonal to the boundary surface. In 
addition, at the corners the boundary data has to satisfy 
compatibility conditions. We now repeat the same anal- 
ysis for the semi-discrete system in order to determine 
appropriate boundary conditions for the computational 
grid. In particular, we examine the application of bound- 
ary conditions at the corner points of the grid. 

We define the following one dimensional scalar prod- 
ucts between vector valued grid functions, 



Ni 



(u, V)h2 



N2 

i2 ^ ufviai , 

1=0 i=0 

(17) 

where ai = {1/2, 1, . . . , 1, 1/2}. The 2D scalar product 
is 



Afi N2 

(u, v)h = hih2 ^ ^ ufjVijaiaj . 
■i=0 j=o 



(18) 



To simplify the notation we introduce = D^^''. If we 
approximate di with the second order centered difference 
operator D^^^ij ~ [ui+ij — Wi_i.j)/(2/ii) in the interior 
(1 < i < TVi - 1, < j < iV2) and with the first order 
one-sided difference operators D^^'^uoj ~ (uij — UQj)/hi, 

D-'^UNij ~ (""A^i.i ~ UNi-i,j) /hi at the x"^ ~ const, 
boundary we have that 



(1), 



(19) 



N2 / Ni Ni \ 

j=0 \ 1=0 i=0 / 



Similarly, if D^^) D^^^ in the interior and D^^) = D^f' 
at the = const, boundary, we have that 



(20) 



If these simple finite difference operators are used to 
approximate the spatial derivatives in, for example, (|15|l . 
the time derivative of the discrete energy 



E = {u,Hu)h = hih2 uf.jHijUijaiaj 



(21) 



gives 
d 



dt 



-E 



+ iu, (dtH + HB + [HBf - d^{HA'))u)h ,(22) 



where we have not assumed energy conservation. 

According to the discrete energy estimate above, to 
control the energy growth due to the boundary term, 
one should give data to the incoming variables in the 
direction ri, orthogonal to the boundary in maximally 
dissipativc form, as shown in Fig. [3 To define the unit 
normal at the corner of the grid we examine the contri- 
bution to the energy estimate due to the corner point 
itself ll^. We see that, for example, at (i,j) = {Ni,N2) 
we have 



^U^N^N2iHA^u)N,N2 + Y^N^N2(HA^u)n,N2 
= ^-J-'^NiN2iH^'^'^)NiN2, 



(23) 



where \h\ = ^h\ + h\ and n = {h2,h\)l\h\ is the unit 
normal at {Ni,N2)- Similar results hold at the other 
corners. In particular, this shows that for uniform grids 
[hi = /i2), data should be given in the 45° direction. 
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FIG. 2: The energy estimate for the semi-discrete initial- 
boundary value problem on domains with corners shows that, 
in order to control the growth due to the boundary term, 
boundary data must be given to the incoming modes with 
respect to the unit normal n. At the corner, the unit normal 
depends on the mesh spacings hi and /12. 



8. Olsson's boundary conditions 

Let's assume that at the boundary there is one incom- 
ing, one outgoing; and one zero speed mode and that 
A = diag(+A+, -A_, 0) with A± > 0. At each grid 
point belonging to the boundary, boundary conditions 
are implemented according to Olsson's prescription t29j. 
Namely, if n = (ni, 712) is the outward pointing imit nor- 
mal, we carry out the following steps: 



1. Compute (W^- 



Id 



Id 



(0;ri)NT 



old 



Q{n)'^Il, where 11 is the discretized right hand side 
and Q{n,) is the orthogonal matrix that diagonalizes 
the boundary matrix Q{n)'^ HA^Q{n) = A. 

2. If the boundary condition at the continuum is 
— -|- g, ovcrwritc the ingoing 



w 



7 



and outgoing modes according to 



new _|_ ^2 ^ old 

1 (+A+;n) 



(-A_;n)^ 



1 



old 



w(-A-;n) 



l + S'2 



Id 



(-A_;n)^ 



old 



l + S'2 
S 

1 + 52 



5*3 
9tg 



and leave the zero speed mode unchanged, 

y* new 

Swii^"'"' +9t.g and that the foUowing Unear com- 
bination of in- and outgoing modes remains un- 



W^i^^'l This will ensure that M^i^w+'"^ 



changed, SW^i+w^'"^ 



W, 



i-X-:n) 



old 



^oid^ Note that unless S = 0, the outgoing 
mode will be modified. When the exact solution 
is known, the boundary data required to reproduce 
it are g = - Sg'^-^~'"\ where 

and (^(~^-'") are ingoing and outgoing characteris- 
tic variables of the exact solution. 

3. The new modified rhs is obtained by multiplying 

(M^i+w+^"\lvi-w-^"\H^ic°w"^r by Q{n). 



9. Consistency at corners 

Although giving data to the incoming variables at the 
corner in the direction n controls the energy growth and 
therefore ensures numerical stability, to achieve consis- 
tency with the boundary conditions used at the two ad- 
jacent sides some extra care is required. Let us assume 
that the normals to the two sides defining the corner are 
n and rh and that A = diag(-|-A+, — A_, 0) with A± > 0, 
i.e., on each side there is one ingoing, one outgoing and 
one zero speed mode. We give data to the incoming vari- 
ables at the sides according to 



^(+A+;n) ^ („) 
'"new y I 

,,,(+A+;ni) _ „(m) 
'"new J ' 



(24) 
(25) 



where, for simplicity, we have assumed no coupling to 
the outgoing fields. At the continuum these two condi- 
tions will be satisfied also at the corner. Let us assume 
that at the corner data is given in the direction p. We 
must translate (|24|) and H25|l in terms of characteristic 
variables in the direction p. If Q{r) denotes the orthog- 
onal matrix defining the characteristic variables in the 
generic direction r, w^"^^ — Q^{r)u, then we find that at 
the corner we must use ly^-'^+'P) = Siu'^'^^'P^ + g^'''\ with 
a non-trivial coupling 



S 



\m)Q{p)],4QTin)Qip)] 



{n ^ m) 



[Q^im)Q{p)]i3[Q^{n)Q{p)]n - (n ^ m] 
and boundary data 

[Q^(m)g(p)]i3ff(")-(n^m) 



[Q^im)Q{p)]i3[QT{n)Qip)]ii - {n ^ m) 



, (26) 



, (27) 



where [Q]ij denotes the ij matrix element of Q. The 
notation (n <-!■ m) indicates that the preceding term is 



repeated with the exchange of the vectors m and n. In 
particular, if g^"^ and g*-™^ vanish, then g^P") also vanishes. 
However, in general, the absence of coupling on the two 
adjacent sides is not consistent with a vanishing S at the 
corner, Eq. 



B. The massless scalar field on a curved 
background 

1. The axially symmetric system 

We now turn to the massless scalar field propagating 
on a curved background {M, g). The equation of motion 
is the second order wave equation 



0, 



(28) 



where V denotes the covariant derivative associated with 
the Lorentz metric g. In terms of the tensor density 

1^" = 



^ the wave equation can be written 
S;.(7''"a,<i>) = 0. 



(29) 



We introduce the auxiliary variables T = dt^ and di = 
and rewrite Eq. (|29|) in first order form. 



dt^ = T, 



(30) 



dtT = ~ (7"9,r + d,{-i''T) + d,{j'^d,)+ (31) 

+ da"T + da'^dj) 

dtd, = d,T. (32) 

The $ component of a sufficiently smooth solution 
(<I>,r, d;) of the first order system satisfies the second 
order wave equation provided that the constraints Ci = 
di — di^ = are satisfied. An attractive feature of this 
particular first order formulation is that the constraint 
variables propagate trivially, namely dtCi = fl3|. In 
particular, this ensures that any solution of (|30II32|) 
which satisfies the constraints initially, will satisfy them 
at later times, even in the presence of boundaries. 

Since $ does not appear in Eqs. H31II32|I . we will drop 
Eq. H3U|) from the system. The constraints are replaced 
by Cij = d[idj] = 0, which also propagate trivially. In- 
terestingly, if Eq. H32|l and the constraints are discretized 
using difference operators satisfying [Di,Dj] = 0, which 
is usually the case, then the time derivative of the dis- 
crete constraint variable Cij — D^idj-^ will also vanish. In 
particular, for initial data such that di = 0, the discrete 
constraints will be identically satisfied during evolution. 

To simplify the problem we assume that the back- 
ground metric is axisymmetric, which implies that there 
exists a spacelike Killing field tp = ijj^^di^ = d^. We al- 
ways use coordinate systems adapted to the Killing field, 
so that the metric components are independent of the 
4> coordinate and, in particular, d^j'^''' = 0. Since we 
are only interested in axisymmetric solutions of the wave 
equation, i.e., solutions which do not depend on 0, the 
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variable can be eliminated from the system. Thus, 
the first order axisymmetric wave equation consists of 
Eq. H^H) and (|32() . where the Latin indices now span only 
two dimensions, and one constraint. 



2. Characteristic speeds 

The characteristic speeds in an arbitrary direction n, 
with |7i| = 1, are given by the eigenvalues of 



A" ^ A'n, 



These 



eigenvalues 
^)/(- 



-27*"/7* 

are 

■7") = 







(33) 

(7*" ± 
and So = 0, 



where a is the lapse function, the shift vector, and 
hij is the induced 3-metric on the t = const, slices in 
the ADM decomposition For the system to be 

hyperbolic it is essential that (7*")^ _^tt^"" = }fin > 
which will be true as long as the t ~ const, hypersurfaces 
are spacelike. We also need s± to be bounded in the 
domain of interest, which will be the case in a cylindrical 
or spherical coordinate system (for r > rg > 0), provided 
that the solution does not depend on the azimuthal 
coordinate 6. 



We assume that in a neighborhood of the outer boundary 
dt is timclike. The integrand of the surface term can be 
written as 

2{Tj*'T+TY^dj)n, = A+u;(+^+'")^-A_?i;(-^-'")^ (38) 
where A± = 7" ± 7*" and 

,„(±A±;») 



w 



Jl ± 7*" 1 7™dj , , 

V2 V2y/1± 7*" 

= 7lrf« (40) 

are the orthonormal characteristic variables of HA"' . To 
simplify the notation wc have introduced the quantities 
7" = y3~7^™7^, Y"" = 7^"/7" and The latter 
satisfies (5ij7^7;^ = 1 and (5y7^7"'" = 0. To express the 
primitive variables in terms of the characteristic variables 
we invert Eqs. Ip^ - lPHjl . 



T 




(-A-;n) 



(41) 



+ 7lu;(°'")(42) 



Equations (|39|I - H42|I will be used in the boundary condi- 
tions. 



3. Symmetrizer, conserved energy, and characteristic 
variables 

One can verify that 

H{t,x)^v{t,x)(^'f ) (34) 

is the most general symmetric matrix that satisfies 
HA"^ = (HA^)^ . When positive definite, which will be 
the case if and only if dt is timelike and 77 > 0, it rep- 
resents the most general symmetrizer of system (|31l I32II . 
If we use a coordinate system adapted to the timelike 
Killing field k = dt, the components of 7^"^ will be time 
independent. In this case the symmetrizer 

H = ( ) (35) 

satisfies Eq. (fT^ and gives rise to a conserved energy. 

The boosting of the black hole will be performed by 
a Lorentz transformation. The time coordinate of the 
boosted frame will no longer be adapted to the timelike 
Killing field. In this case the time derivative of 

E = [ {--f'^T^ + f^d,dj)d^x (36) 

is given by 

^ 2 f {Tj''T + TY^d,)n,da (37) 
Jon 

+ [ {Tdt-/"T + 2Tdti''d, + dAi'^dj) d^x 
Jn 



4- Discretization 

Even when there is no conserved energy, it may be 
desirable to discretize the right hand side of 1)31(1 and 
(|32(l in a manner that satisfies the optimal estimate H22|l . 
such as (|16f) or other alternatives, where the symmetrizer 
H is given by 

The discretization of the wave equation according to 
(Unj leads to 

dtT = -(^^*'D,T + D,{f'T) + ^D,{Y^d/)+ 

+ ^f'D,d, + ia,7'^d, + dt^"T + dti'^d,^ /7", 

^td^ - iAT+i(V).fe^,(7'^r)-i(V).fe9,7'^r, 

where •^7"^ denotes the inverse of 7*-' . 

Alternatively, one can simply replace the partial 
derivative dj in Eq. 1(31(1 and 1(3 2() with the finite differ- 
ence operator Di satisfying l(19() and ((20(1 and obtain the 
semi-discrete system 

dtT = -{j'W,T + D,{f'T) + D,{f^d,)+ (43) 

+dtj''T + dt-/''dj) 

dtd, = AT. (44) 

which also satisfies the estimate ((22(1 . It is this discretiza- 
tion that will be used throughout this work, even in the 
boosted black hole case (9(7^" 7^ 0), where the energy 
(136(1 is not conserved. We analyze the discretization at 
the axis of symmetry in the next sections. 
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IV. MINKOWSKI BACKGROUND 



3. The boundary conditions 



The energy method for constructing stable finite dif- 
ference schemes has, until recently, received little at- 
tention in numerical relativity. Thus wc first present 
the wave equation in axisymmetric Minkowski space to 
demonstrate the method, before moving to the more com- 
plicated black hole configurations. In this section we 
give energy preserving discretizations for cylindrical and 
spherical coordinates. In particular, we will show how 
to discretize the system on the axis of symmetry in an 
energy conserving way. The next section examines dis- 
cretizations for a Schwarzschild black hole in Kerr-Schild 
spherical coordinates. 



Since in this coordinate system dt is a Killing field, the 
energy H36|l is conserved. The time derivative of 



E 



(T^ + P^ + Z^) pdpdz (51) 



gives only boundary terms which can be controlled by 
giving appropriate boundary data 



dt 



-E ^ 2 



T{t, Pmax, z)R{t, Pmax, ^)Pmaxrf2 (52) 



^2 / [T{t,p,z)zit,p,z)]iziz::pdp. 



A. Cylindrical coordinates 



4- Energy conserving discretization 



1. The system 

In a Minkowski background in cylindrical coordinates, 
{t, p, z, <^}, the second order axisymmetric wave equation 
has the form 



We consider the first order formulation 



(45) 



The discretization of the right hand side of Eq. 1)46(1 
at the p = axis deserves special attention. As a conse- 
quence of the regularity conditions we have that 



hm -dnipP) 
p_0+ p '"'^ ' 



2dpP\ 



p=0 ■ 



(53) 



and therefore no infinities appear on the right hand side. 

This suggests considering the semi-discrete approxima- 
tion ii 



dtT = ^dp{pP) + d,Z, 
dtP = dpT, 

dtZ = a,r. 



(46) 

(47) 
(48) 



where T ~ dt^, P ~ dp^ and Z ~ dz^ arc functions of 

{t,p,z) e [0, T] X [0, pinax] X [Zmin, Zmax]- 



2. Regularity conditions at the axis p — 

Smoothness at the p = axis requires that the odd p- 
derivatives of the scalar field vanish on the axis, namely 
a2"-i$(i,p, z)|p=o = for n 1,2, . . .. This implies 
that the following conditions for the auxiliary variables 
T, P, and Z, have to hold during evolution 



PI 



p=0 



9„'"P|p=o = for = 1,2,... (49) 



dl"-'T\p=o 



dT~'Z\p^, 



for n = 1,2,.. (50) 



If the initial data satisfies H49|l and H50() . and the pre- 
scription P{t, 0, z) = is used as a boundary condition 
at yO = 0, then the above conditions will hold at later 
times. 



2D[p'^P 



Oj 



D'^'-^Z, 



Oji 



\ l-D(p\pP),j+D('^Z, 
= D'-p'^T,j, i > 1 



i = 
i > 1 



(54) 

(55) 
(56) 



where pi = lAp and Zj = Zmin+jAz, with NpAp — pmax 
and NzAz = Zmax — -Zmin- Thc difference operators 
D'-p^ and are second order accurate centered dif- 

ference operators where their computation does not in- 
volve points which do not belong to the grid, and are first 
order accurate one sided difference operators otherwise. 
The regularity condition, Poj = for j = 0, . . . , A^z, is 
enforced for all t, and Eq. (|49|) ensures that D^^'' Pqj is, 
in fact, a second order approximation. A solution of 1(541 
1551 156(1 conserves the discrete energy 



1 



E (^5 + + 4) P^^^^P (57) 

(7jAz , 



1=1 

.2 



which is consistent with thc continuum expression 1(51(1 . 
More precisely, using the fact that dtToj = -^Pij + 

D^^^Zoj and the basic properties of the finite difference 
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operators, one ean sec that the following estimate 
d 

-E = 2Y,Tn,jPnPn,,<Tj^z (58) 

N, 
i=\ 

holds, consistently with the continuum limi t (I52|l . 

As it is pointed out in Section 12.7 of [23, one or- 
der less accuracy at the boundary is allowed, in the sense 
that it does not affect the overall accuracy of the scheme, 
provided that the physical boundary conditions are ap- 
proximated to the same order as the differential operators 
at the inner points. 

5. Discrete boundary conditions 

By inspecting the boundary terms of the discrete en- 
ergy estimate (|58|) . we can readily see how the bound- 
ary data should be given at each boundary grid point 
(j = 0, j = Nz, and i = Np). In the case of a uniform 
grid (Ap = Az), in order to control the energy growth 
boundary data should be given in maximally dissipative 
form in the directions shown in Fig. |21 



on the axis, indicates how to specify boundary data at 
the corner grid points that lie on the axis. 

B. Spherical coordinates 

In this section we discretize the wave equation in 
Minkowski space with spherical coordinates {t,r,6,(j)}. 

1. The system 

The second order axisymmetric wave equation on a flat 
background in spherical coordinates 

= }-dr{r'^dr^) + —^deisinOde^), (59) 
sm ti 

is written in first order form as 

dtT = lo,(r2i?) + — ^ae(sin0e) (60) 

dtR = drT (61) 
dtQ = deT, (62) 

where T — dt^, R — dr^ and 8 = are functions of 
(t,r,0) e [0,T] X [r,„i„,w] x [0,^]. 

2. Regularity conditions on the axis 9 — and 6 — -k 

Smoothness requires that the odd 6'-derivatives of the 
scalar field vanish on the axis of symmetry, namely at 
6* = and 6 = n. This implies that 

e\g=o,n = d^g''e\g=o,. = n = 1,2, . . . (63) 
-^T|e=o,. - -ii?|e=o,. - n= I, 2, .(.64) 

As in the cylindrical case, it is possible to show that if 
the initial data satisfy H63|) and (|64|l and the boundary 
conditions Q\e=o.Tr — are used during evolution, then 
the above regularity conditions will continue to hold. 



3. Boundary conditions 



FIG. 3: This figure shows how the unit normal at the bound- 
ary grid points should be chosen. We note that at the corners 
which lie on the axis of symmetry we must apply both the 
regularity and the boundary conditions. 



The presence of lower order terms in ()57f) , in addition 
to ensuring that the discrete energy is positive definite 



Since we are interested in a domain of the form ~ 
{{r,e) e R^lrmin < r < rmax,0 <9<Tr}, where r^in > 
0, the characteristic speeds are bounded by max{l, r~^l^.^}. 
The conserved energy is 



E 



92 

2 



/o V r 
and its time derivative is given by 

d ^ " ' r^RT\ 



svaOdOdr, (65) 



dt 
and r 

to the incoming modes. 



sin Ode. 



(66) 



boundary data must be given 
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4-. Energy conserving discretization We discretize the right hand side of (|60H62(I as 

As a consequence of the regularity conditions on the 
axis of symmetry, we have that 

Um -^deismee) = 2dee\e=m^, m = 0, 1. (67) 

B^miT Sm u 



dtT, 



dfRij 



1 ^D(-\r^R),,+ 



-D^'' (sin 09),, J 



(68) 

(69) 
(70) 



where = rmin+i Ar and 9j = jA6, with NrAr = rmax" 
Tmin and NeA9 = it. The condition O^o = QiNo = on 
the axis is enforced at all times. The following discrete 
energy 



r>2 



Nr Ne-1 

^ = EE 

4=0 J = l 

+ lj2iT"o + Rlo)rla.sinA9^Ar 

i=Q 

1 A9 
+ 2 E i^^No + ^k) ^''^^ sin A9—Ar, 

i=0 

is conserved by the semi-discrete system. Its time deriva- 
tive 



(72) 



{TiNe R 



\i=N, 



smA6' , 



gives only boundary terms consistently with the contin- 
uum estimate (jBBJ. 



5. Discrete boundary conditions 

The choice of unit normal at the boundary grid points 
{i = and i = Nr) is illustrated in Fig. 21 



V. FIXED BLACK HOLE BACKGROUND 

This section is a generalization of the results of the 
previous section to the case of a static black hole back- 
ground. The background metric is Schwarzschild in 




FIG. 4: According to the discrete energy estimate, boundary 
data should be given to the incoming variable in the direction 
indicated in the figure. 



Kcrr-Schild coordinates [43 . The Cartesian components 
of the background metric can be written as 

2M , , 

J. f 

where 77^,^ = diag(— 1, -1-1, -1-1, -1-1), r'^ = + y'^ + and 
~ {l,x/r). In these coordinates the determinant of 
the four- metric is g = — 1. 

The tensor density components 7'"', which are needed 
to write down the 3D wave equation in first order form, 
are given by 



r 



T]' 



2M 



(74) 
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where = ■q'^'U = (-l,x/r). 

As we do not wish to consider cyUndrical excision in 
this paper, we analyze here only the spherical coordinate 
case. 



that 



n = 1,2, . 



(79) 



T\e=o,. = d'/'~'R\e=o,n = n = 1, 2, .(.80) 



A. Spherical coordinates 

1. The system 



3. The boundary conditions 



In spherical Kerr-Schild coordinates the components 
of j^" on a Schwarzschild background are 



sin t 



diag 



-i,+i, 



2M 
r 
1 



(75) 



1 



= (-1, +1,0,0). 
The first order axisymmetric wave equation is 

-ae (sin 616) 



2 M 9 M 1 

—drT + —dr{rT) + —dr{rr-R) (76) 



1 



dtR 



rr^ sm ( 

= drT 



(77) 
(78) 



where r± = r ± 2A/, T = (9t$, i? = 9^$ and 6 = ^e*. 
In the region of interest, n = {(r,6') e R2|2M < r < 
''maxjO < < tt}, with A/ > 0, the characteristic speeds 
are bounded. 



A symmetrizer which gives rise to a conserved energy 
is given by 



H = diagjrr^ sm9, rr sm9, sin 9}, 



(81) 



which is positive definite for < 6* < tt and r > 2M. 
Inside the event horizon, r < 2M, the vector field dt 
becomes spacelike and the system in only strongly hy- 
perbolic. 

The time derivative of the energy 



E 



{rr+T^ + rr-R^ + 6^) sin9d9dr, (82) 

l2M Jo 

gives only boundary terms 



dt 



-E 



(2MrT^ + rr- TR) ' 7!" sin 

V / I r—2M 



(83) 



In addition to the regularity condition O = at the 
axis, the problem requires boundary data at r = fmax. 



2. 



gularity conditions on the axis 9 = and 9 = n 



Smoothness requires that the odd ^-derivatives of the 
scalar field vanish on the axis of symmetry. This implies 



4-. Energy conserving discretization 
We discretize the right hand side of H76l - I78|l as 





( 2M^0 






dtRij 









2M 
2M 



1 

ViT ■ 
1 



-D''''\rr-R)ij 
-D^'''>{rr~R)ij 



' 1 



(9) 



D)^\sin9Q), 



0, A^e 

1 Na 



(84) 

(85) 
(86) 



where n = 2M + iAr and 9j = jM, with iV^Ar = -2^ 

rmax-2M and NeM = tt. The following discrete energy, +3 ^^o + nr, R,o) a^ smA9—Ar 

4=0 

1 A9 

Ne-1 +2 H i^i^tTMg + nr-RiNe) sin A6I— Ar, 

^ = E E ("-^"-tT^j + nr-Rl + el) sm9,a,A9Ar 

i=0 j=l 
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is conserved. Its time derivative is given by 

= 2 ^ (2A/r.7;5+nrrr,,i?,,)|;:^'-sin0,A0 
i=i 

+ {2MnTlo + r,r^T,oR,o) sin 

+ (2Mr,T2^^ +r,r-r,jv«i?.jv«)|;:J'-sinA0^. 

We point out that, since 8^0 = QiNg = 0, the discrete 
energy is positive definite on the axis. However, because 
Tq = 0, it does not control the growth of Roj on the 
event horizon. Numerical experiments indicate that this 
does not cause any problems. Moreover, experiments do 
not suggest that placing r^in within the horizon, where 
the equations are only strongly hyperbolic, leads to an 
unstable scheme. 



5. Discrete boundary conditions 

Data should be given to the incoming modes at r = 
Tmax as in Fig.0] Unlike the Minkowski case, no bound- 
ary conditions should be given when the inner boundary, 
r = rniin, is at or within the event horizon. 

VI. BOOSTED BLACK HOLE BACKGROUND 

Finally, we consider the case in which the scalar field 
propagates on a boosted black hole background. To solve 
this problem we introduce two coordinate patches: one 
patch fixed to the outer boundaries and one patch co- 
moving with the black hole, and fixed to the inner exci- 
sion boundary. We choose cylindrical coordinates for the 
first coordinate patch, boosted with respect to the black 
hole such that the hole moves with velocity /? along the 
symmetry axis in these coordinates. Spherical coordi- 
nates are used on the second patch. These coordinates 
are chosen by fixing the event horizon at a constant co- 
ordinate value (r' = 2M), and requiring that all data in 
both coordinate systems are simultaneous. By adapting 
these coordinates to the black hole horizon, we may ex- 
cise the spherical grid at r' = 2M for all values of the 
boost parameter. In this section we first write down the 
the components of the 4-mctric in a boosted Cartesian 
coordinate system, then discuss the two coordinate sys- 
tems in some detail. 

We recall that in a Cartesian coordinate system 
{t, X, y, z}, with respect to which the black hole is at rest, 
the metric components have the form given in H73|l . Un- 
der a Lorentz boost, i.e., in the new coordinates 

t = 7(f - /3z) (87) 

X = X 

y = y 

z = j{z-f3t), 



where 7 = (1 — /3^) the components of the Kerr- 
Schild metric become 

2M 

?7^p diag{-l,+l,+l,+l} , 
^^L = {f,x,y,z)/r, 

where f — ^{r + Pz), z — 7(2 -f- /3r), z = 7(2 -1- and 
r"^ = + z^ =x^ ^ ^\z -f /Jf)^. At time t the 

singularity is located at (a;, y, z) = (0, 0, — /3i). 



A. Boosted cylindrical coordinates 



z 




2 4 6 

P 



FIG. 5: The regions delimited by the curves are regions in 
which the system in not symmetrizable hyperbolic, but only 
strongly hyperbolic for /3 = 0,-1/4,-1/2,-3/4 (the black 
hole is located at (p, z) — (0, 0) and moves in the +z direc- 
tion) . As the boost parameter /? increases in magnitude there 
is a larger part of the cylindrical domain in which the system 
in only strongly hyperbolic. 

We now choose cylindrical coordinates {i, p, z, 0}, with 
pcos0 = X and psin^ = y, giving 

'qf.y = diag{-l, +1, +1, +/?^} , 
(■ii = {r,P,z,0)/r 
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and 



P V 



r 



diag{-l,+l,+l,+— }, 
P 

{-f,p,z,0) /r, 



where z = 7(2 + /3f) and = + 7^(2 + Un- 
fortunately, in these coordinates the wave equation has 
a rather unpleasant form: the components of 7'^'' have a 
non-trivial dependence on the three coordinates p and z, 
and especially, i. 

The analytic expressions for the time derivatives of the 
7*'^ components are needed. Using the fact that div = 
P^z/r, dtf = I3^z/r and diz = P^f/r we get 



da'' 



2M/37pf(32f - 2rz)/r^, 
2M/37p2(rz-3zf)/r^ 
2Ml3jp {rz^ + rf"^ - 3zfz) /r^. 



(89) 



In this coordinate system our first order formulation 
has no conserved energy {d^^'^ 7^ 0). The region in 
which the system is symmetrizable hyperbolic is deter- 
mined by the set of points in which dt is timelike, 



1 



2MP 



> 0. 



(90) 



Fig. El shows the regions of lack of symmetric hyperbol- 
icity for different values of the boost parameter. 

On the axis of symmetry (p = 0) the equations need 
to be expressed in a form which avoids "0/0". This can 
be done by taking the limit p in the equations. It is 
convenient to introduce the quantities 



p 



p^ 



r = ^, 

p 



vPP 



^IP^ — 



P^ 



which have a finite limit for p — > (since the singularity 
is excised we can assume that r > ro > 0). The right 
hand side of at p = becomes 

dtf = (7*^9jT + 27*'^T + 9j(f^"T) + 27'''^apP-f 

+27^^^"^ + d-,{TZ) + ditf + d-a^'^Z] /(-7") . 



B. Co- moving spherical coordinate system 

We introduce a spherical coordinate system 
{<', r', & , (/)'} which is related to the unboosted Cartesian 
coordinates {i, x, y, z}, the coordinates in which the 
black hole is at rest, by 



t' = t = 



t = lit-f3z) 



(91) 



y''x^ + + z'^ 



tan 



We emphasize that (|91|) is not a Lorentz transformation. 
The coordinates are adapted to the event horizon in the 
sense that its location (r' = 2M) is time independent, 
and setting t' = i maintains simultaneity in the two co- 
ordinate systems. The metric components and the com- 
ponents of "fP " are more conveniently written in matrix 
form 



(1 



i (ML 
/3 cos 61') (1 



-/3cos0'(l-^)) i(r'- 
/3cos6i' + ^(1 + /3cos6i')) /3sin6'' ((/ - 



2M)f3sint 
2A/)/3cos6 



2M) 



e'ir' - 2M)) 



r'2 sin^ e' / 



/ -ir'sin6l' (r' + 2M-i'^{l +Pcos9')^) r' sinO' {2 M - (3 cos 9' {r' ~ 2M)) f3r' 



' sin^ 0' 



ir'(r' - 2M)sin6'' 





sin 9 



\ 





inW / 



7 sin 



r 



(92) 



We have symmetric hyperbolicity for r' > 2M and < 
0' < n. 



On the axis of symmetry {0' = or 0' = tt) the equa- 
tions need to be expressed in a form that avoids "0/0". 
This can be done by taking the limit 0' 0q, where 
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^0 = 0, TT in the equations. If we introduce the quantities 

„t'i 















sin 6' ' 


sin 9' 








sin 6*'' 




sin 6"' 



sin^ 0' 



which have a finite hmit for sin0' 0, then the right 
hand side of H31|) on axis becomes 

dt'T' = (f '■'a.'T' + a.'(f' '^'r')±2f' ^'r' (93) 

+dr>{r'''' R') + 2^''0'dB.Q') /(-f *'), 

where the components of 7^* are understood to be eval- 
uated at 9' = 0, TT. 

C. Overlapping grids 

As mentioned in Section^ the method of overlapping 
grids gives a natural method for solving finite difference 
problems on multiple domains. For the boosted black 
hole, we use the cylindrical grid as our base grid and 
introduce the spherical grid adapted to the inner bound- 
ary (event horizon). The two grids overlap as shown in 
Fig. El 

The spherical grid requires boundary data at r = r,nax, 
which does not constitute a physical boundary in this 
problem. Here the data is computed by interpolating the 
values of the field from the cylindrical grid. Similarly, the 
values of the fields at the grid points of the cylindrical 
grid near the excision region, which lack a neighboring 
point in a coordinate direction (these points are marked 
with a square in Fig. O, are also updated via interpo- 
lation. In this work we used second order Lagrangian 
interpolation, which, for a scalar quantity, is given by 

hnt{xi + a Ax, yj + bAy) = (1 - a)(l - 6)/„- 
+ (1 - a)fo/,;j+i + a(l - 6)/i+i,j + a6/j+ij+i 

where < a,5 < 1. Higher-order interpolation stencils 
may also be used, though for the cases examined here, 
improvements in the resulting solutions are slight, result- 
ing in no increase in the order of accuracy. 

The boosted cylindrical and co-moving spherical coor- 
dinate systems are related by 



t = t' 
p ^ r' sin 9' 
z = ^~^r' cos 9' - (it' 



(94) 



The evolved fields are not scalar quantities, but com- 
ponents of a 1-form. So, in addition to the coordinate 
transformation between the two coordinate systems, the 
communication of the values of the fields requires the use 
of the transformation law of 1-forms. In this case we have 



T = r + J (3 cos 9' R' --f/3- e' (96) 

r 

p = sin^?'i?' + H^e' 
sin 9' 

Z = 7cos0'i?'-7^e', 



and 



T' = T~(3Z 

, - cos 9' - 
R' ^ sin9'P+ Z 



e' = r' cos9'P-r'^-^Z, 

7 



(97) 



where (T,P,Z) and {T',R',Q') are the fields on the 
cylindrical and spherical grids, respectively. 



1. Discretization on the axis 

The discretization of the system in boosted cylindrical 
coordinates in the interior and at the outer boundary is 
done according to Eqs. (|43|I - H44|) where the components 
7'^'' are given in H88|) . On the axis of symmetry we use 

Q-f = (ji^-D'^'-^f + 2fPf + D(^) (7*""T') + 2YPDfp+ 

2Y"z + D^^^^^'z) + an"T + da'^z) /(-f*). 



Similarly, the discretization of the system in co-moving 
spherical coordinates in the interior and on the event 
horizon is done according to Eqs. (|43|I - H44() . where the 
components of 7^^ are given in (|92|l . On the axis of 
symmetry (9 — and 9 = tt) we use 

dt'T' = (^f''''D^'''^T' + D^'''\f''''T')±2^'''''T'+ 

d^-'HY'^'r') + 27«''''z?i'''e') /(-7*'*'). 

2. Boundary conditions 



and the inverse transformation 
t' = t 

r' = V/o' + 7'(^ + /3tT 



9' 



tan 



-1 



l{z + (3t) 



Boundary conditions in maximally dissipative form are 
given at the outer boundary of the cylindrical grid in the 
(95) directions indicated in Fig. El In the boosted case, in- 
stead of overwriting the right hand side at the boundary 
according to Olsson's prescription, we overwrite the so- 
lution itself. The reason for doing so, is that it avoids the 
tedious task of computing time derivatives of the bound- 
ary data. 
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It is interesting to notice that the outer boundary of 
the cyhndrical grid could become at some points purely 
inflow (s± > 0) for very large values of (3. (We exclude 
the case in which the black hole is outside the outer 
boundary.) At and near these inflow boundary points 
the system is only strongly hyperbolic and the energy 
method fails to give the correct boundary conditions. As 
it is pointed out in |34j . applying maximally dissipative 
boundary conditions to strongly hyperbolic systems can 
lead to an ill posed IBVP. Where the boundary is purely 
inflow we give data to the two incoming fields. Our nu- 
merical experiments (Sec. I Villi indicate that the scheme 
is stable. 



have a magnitude of 



1 



(99) 



which is greater than 1 for /? 7^ 0. For example, for 
/3 = 0.9 the characteristic speeds in the spherical grid 
can be as large as v'l^ w 4.359. In this case a Courant 
factor larger than « 0.46 is likely to lead to numerical 
instability. 



VII. NUMERICAL EXPERIMENTS 



3. Artificial dissipation 

Whereas the single grid schemes do not require any 
artificial dissipation, it is known that overlapping grids 
require explicit dissipation for stability [Siij . To the right 
hand side of the discretized system a term of the form |33 

Q,u = -a (h\{D'-l^D'-l^f + /i3(i^(2)^(2))2^ ^ (98) 

is added. This dissipative operator is modified near the 
outer and inner boundary, as was done in |l3l |. Near 
and on the axis of symmetry dissipation is computed ex- 
ploiting the regularity conditions of the fields. As this 
dissipation has a five-point stencil, wc find that the long- 
term behavior of the code is improved in some cases by 
interpolating two points at all inter-grid boundaries. 



4- Choice of Courant factor 

The fully discrete system is obtained by integrating the 
semi-discrete system with third or fourth order Runge- 
Kutta. Whenever explicit finite difference schemes are 
used to approximate hyperbolic problems, the ratio be- 
tween the time step size k and the mesh size h = min{/ii}, 
the Courant factor, cannot be greater than a certain 
value This Courant limit is inversely proportional 

to the characteristic speeds of the system. We estimate 
allowable values for the Courant factor by examining 
the 2D wave equation in first order form, dtu^ = d^Ui, 
dfUj ~ diUQ. Assuming second order, centered differenc- 
ing for the spatial derivatives, we plot the Courant limits 
for third and fourth order Runge-Kutta as a function of 
the artificial dissipation parameter in Fig. |3 

The characteristic speeds in the cylindrical grid arc 
bounded by 1 in magnitude. Looking at Fig. [3 this would 
suggest that one could use a Courant factor as large as 2.0 
(using fourth order Runge-Kutta and ignoring the fact 
that this is a variable coefficient problem with lower order 
terms and with boundaries). However, in the spherical 
grid the characteristic speeds along the axis of symmetry 



To our knowledge there are no stability proofs for two 
dimensional hyperbolic problems with overlapping grids. 
To check the convergence of our code we must rely on 
numerical experimentation. 

Let u{t, x) be the exact solution of the continuum prob- 
lem and vfj the solution of the fully discrete approxima- 
tion. If, for any t, as nk t, 

el = (/ii/ia^I - u{nk,x,,)f fl^ = Om+0{k^) 

a,s k,h ^ 0, the difference scheme is said to be convergent 
of order {p, q) . This implies that the overall order of 
convergence of the scheme, assuming k/h ^ const., is 

e" 

g = lim log2 ^ min{p, q} (100) 

as nk — > t, where t is some fixed time. To use this equa- 
tion we must know an exact solution of the continuum 
problem. 

Exact solutions for the scalar wave equation in 
Minkowski space are well known. To test our overlap- 
ping grid system we use spherical waves |45j given by 

'S> = Y.M^)Pdcose)e-''^\ (101) 

e 

where Pe{cos9) are the Legendre polynomials and we 
choose fi{r) to be the Hankel functions, which asymp- 
totically represent in- and out-going waves. We tested 
the long-term behavior of our code by evolving an in- 
going spherical wave exact solution, and computing the 
norm of the error. These results are shown in Fig. |S1 

When an exact solution is not available, which is often 
the case, the following standard technique of numerical 
analysis can be useful. Let w be an arbitrary function and 
let us rewrite the partial differential equation as L{u) ~ 
0. If w is inserted into the equation, in general, it will 
produce a non vanishing right hand side, 

L{w)=F. (102) 

Clearly, the modified equation L{u) = L{u) — F = has 
w as an exact solution and the convergence of the code 
can be tested using Eq. (|100|l . 
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We chose w{t, r, 6) = sm{t + r) cos{n6), where {t, r, 9} 
are the spherical coordinates of the rest frame and n is 
an integer. This is an exact solution of 

V^,Vw- F ^0. (103) 

where F is given by 

cos TlO 

F = — {2r cos{t + r) - n'^ sm{t + r)) (104) 



Both w and F are scalar quantities. The evolution equa- 
tion 131|) is modified according to 

dtT = - (7"a,r + d.iY'T) + d,{Y'd,) (i05) 
+da''T + da'^dj - V^F) /7", 

where g = det{g^,y). On the axis of symmetry we use the 
limits lime-»o^i^ = n and lim^^^ sini^ = (-l)"+in. 

^ suit? u^-n siYi8 ^ ' 

The results of our convergence tests for different values of 
the boost parameter and n = 1 are summarized in Fig.|51 
Movies are also available [46l |. 

VIII. CONCLUSION 

Systems with moving boundaries arise in a variety of 
situations, and their solution typically involves introduc- 
ing coordinates adapted to the boundaries. These may be 
either a single, global coordinate system, such as those 
used in binary black hole evolutions [i^. Ell that keep 
the black holes and the outer boundary at fixed coor- 
dinate positions, or, as advocated here, multiple coor- 
dinate patches. Whichever approach is adopted, fixing 
coordinates to the boundaries allows one to unambigu- 
ously specify proper boundary conditions, as required for 
well-posed problems. Moreover, boundaries at fixed grid 
coordinates eliminate the need for extrapolated data at 
points that emerge from a moving boundary. 

In our model problem we have evolved an axisymmet- 
ric scalar field on a boosted Schwarzschild background. 
We used a cylindrical coordinate patch with respect to 
which the outer boundary is fixed, and we introduced an 
overlapping spherical coordinate patch co-moving with 
the hole. At any given time in the co-moving coordinate 
system the location of the inner boundary (the horizon) 
corresponds to a r = const, surface. This surface can be 
represented exactly on the numerical grid and allows one 
to smoothly excise a large volume of spacetime, much 
larger than that permitted by cubical excision. Our nu- 
merical implementation made use of overlapping grids, 
where different but equivalent problems are solved on 
separate grids. To communicate data between grids we 
used interpolation on all fields. 

The discrete version of the energy method, based on 
differencing operators that satisfy the summation by 
parts property, has demonstrated to be particularly effec- 
tive for the construction of a stable discretization scheme 



on the axis of symmetry, and for the identification of dis- 
crete boundary conditions. The stability proofs, which 
hold on individual grids, cannot be immediately extended 
to the overlapping grid scheme, due to the interpolation 
of data from one grid to another. We note, however, that 
it may be possible to define orthogonal projection opera- 
tors for the interpolation that could allow to analytically 
demonstrate numerical stability [j^ . This is a question 
of active interest. Nevertheless, our numerical tests in- 
dicate that our scheme is convergent, even for very high 
values of the boost parameter, and does not suffer from 
long term power-law or exponential growth in the error. 

The model problem that we presented in this work 
is primarily intended as a proof of concept, and several 
avenues of research remain to be explored. Foremost 
might be the addition of the Einstein equations to the 
system for a dynamic black-hole spacetime, where the 
locations of the black hole singularity and event horizon 
are not a priori known. Depending on the dimensionality 
of the problem, one or more coordinate patches adapted 
to the inner boundary would have to be generated dur- 
ing evolution, along with the relationship between the 
various coordinate systems. By monitoring the charac- 
teristic speeds on the excision boundary (with respect to 
the coordinate system adapted to that boundary), one 
can guarantee its purely outflow properties, an essential 
requirement of excision. 

Although alternative numerical approaches may be 
possible, the overlapping grid method has struck us for 
its strength and its simplicity. Owing to its flexibility 
in representing smooth, time dependent boundaries, we 
believe that this technique, or a similar one, will play a 
significant role in the solution of the binary black hole 
problem. 
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Appendix 

To gain some insight into the Hmitations of cubical 
excision and the consequent need for a smooth excision 
boundary, we consider the analytic Schwarzschild and 
Kerr solutions [20I l2ll | . We employ the commonly used 
Cartesian Kerr-Schild coordinates, which are smooth 
across the horizon, and write the metric as 

where rj^v is the Minkowski metric, if is a scalar, 

+ cos^ 9 ' 

and ifj, is a null vector, 

^ f rx ^ ay ry — ax z\ 
^ \ ^ r'^ + a'^ ^ r-^ + a"^ ^ r J 

The parameter A/ represents the mass of the black hole 
and J = aM is the total angular momentum. The 
spheroidal coordinates r and 9 are given by 

2; 

cos 9 = - 
r 



where — x"^ + y'^ + . 

We center a cube of side length L = 2b on the black 
hole, x^ G [— b]. In order to excise this region from the 
computational domain, we must ensure that its bound- 
ary is purely outflow, i.e., that no information can en- 
ter the computational domain. To determine the al- 
lowed values of b, we calculate the characteristic speeds 
on each face of the cube and check that the inequal- 
ity s± — 0' where n is the outward unit normal to the 
boundary, is satisfied. The Schwarzschild solution is ob- 
tained by setting a — 0, and the calculation gives 
< 5 < 2V3/9M w 0.385A'f . The calculations for Kerr 
(a 7^ 0) are more involved, and we present our numeri- 
cally generated results in Fig. ^| We find that because 
of the ring singularity {p = a, z = 0), in addition to a 
maximum size for the excision cube, there is also a mini- 
mum size. In addition, we notice that no cubical excision 
is possible for a > 0.0851Af . This is a severe constraint 
on the spin parameter, and precludes cubical excision for 
interesting values of spin. We note, however, that this 
limitation is coordinate dependent and that it might be 
possible to choose coordinates in which cubical excision 
may be done for higher values of a. 
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FIG. 6: The first figure shows the overlapping grids for the 
axisymmetric wave equation on a boosted black hole back- 
ground. The spherical grid is used to excise the black hole 
fi-om the computational domain. The dashed line represents 
a r' = conts. curve on the spherical patch. This is where 
the cylindrical coordinate system terminates. The arrows on 
the outer boundary indicate how the unit normal is chosen at 
each boundary grid point where boundary data must be spec- 
ified. The circles and the square are used to mark points of 
the spherical and cylindrical grid respectively which are up- 
dated through interpolation. As shown in the second figure, 
the value of the fields on the grid point P of the spherical grid 
are computed by interpolating the values of the fields on the 
four neighboring points. 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 



FIG. 7: The Courant limits for the 2D first order wave equa- 
tion dtUo = d'ui, dtUi = diUo for fourtli order (4RK) and 
third order Runge-Kutta (3RK). The calculation assumes 
no boundaries and second order centered difference opera- 
tors for the approximation of the spatial derivatives. The 
value of a represents the amount of artificial dissipation, 
— (T^^^-i^ h^{D^^ D^^Yu, added to the rhs of the equations. 
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FIG. 8: These long term runs done at a fairly course reso- 
lution for j3 — M = suggest that the interpolation between 
the overlapping grids does not introduce any power law or 
exponential growth. Here the exact solution is given by the 
the real part of the i = mode with uj — 2, with the in-going 
Hankel function for the radial variable in Eq. llUlll . The 
dissipation parameter is set to cr = 0.02. The ranges of the 
dependent variables are < p < 10, — 10 < ^ < 10, 1 < r < 5, 
and < < n. The coarsest resolutions for the cylindrical 
and spherical grids are 90 x 170 and 50 x 68, respectively. 



22 



4.0^ 

3.5 - 
3.0 - 



cylindrical 
spherical 



10 12 



I 

14 



cylindrical 
spherical 



O) 2.0 ■ 



1.5 ■ 



1.0 ■ 



0.5 



O.OL 



1 2 3 4 5 6 
t 



10 11 



4.0^ 
3.5 - 
3.0 - 



cylindrical 
spherical 



FIG. 9: In these figures we plot the convergence factor 
Q{t) = log2(eh/eh/2) a-s a function oi t = nk for a sufficiently 
small value of the spacing h. The convergence test was carried 
out with the modified system Hl()5^ for boosting parameters 
/3 — —0.5, —0.75, and —0.95, from top to bottom, suggesting 
that the equations are correctly implemented and that the 
overall scheme is second order accurate. The resolutions used 
in the cylindrical and spherical grid are 256 x 512, 128 x 384 
and 512 x 1024, 256 x 768. The domain extends from -lOAf 
to +10M in the z direction and up to +10M in the p direc- 
tion. The spherical patch covers the region 2M < r' < 3M. 
We used a Courant factor of 1.15, 0.75, and 0.32, and a dissi- 
pation parameter of cr = 0.02. The evolution is stopped just 
before the spherical grid touches the outer boundary of the 
cylindrical grid. In the /? — —0.95 the system in only strongly 
hyperbolic at the bottom right corner of the cylindrical grid. 
We found that, to achieve convergence, we must give data to 
all fields at this point. 
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FIG. 10: This figure indicates the limitations of cubical exci- 
sion in the Kerr spacetime in rectangular Kerr-Schild coordi- 
nates. We assume that the excision cube is centered on the 
hole, and that the faces of the cube are at ±&. (See descrip- 
tion in text.) Values of 6 for which an inner boundary has no 
incoming modes, and thus a candidate for an excision bound- 
ary, are indicated by the shaded region in the figure. The 
structure of the Kerr spacetime results in both maximum and 
minimum limits to the size of the excision cube. We find 
that, in this particular coordinate system, cubical excision 
for Kerr is well-defined only for very small spin parameters, 
a < 0.0851M, where a = J/M. For values of b below the 
dashed line, the excision cube intersects the ring singularity. 



